
** Spatial approach, 3: Define Neighborhood ** 

		**********------------------------------**********
	* Contents: 
		* Build weight matrices for a given neighborhood definition
		* Beam corresponds to "Strips" in final paper  
		* Input: cell-level data based on spatial_2_collapse_s 
		* Output: weight matrixes in .dta format 

capture set maxvar 30000, permanently // required to run file.  
	
capture log close 
clear all
set more off 

cd "R:\WSV2\TBu_AKe\Spatial_NEW"

capture mkdir Data
global store "R:\WSV2\TBu_AKe\Spatial_NEW\Data"

log using spatial_3_wmat_beam, replace 


cd "R:\WSV2\TBu_AKe\Spatial_NEW"


use "C:\Users\hy65byfe\Desktop\smerge_0712\spatial_2_collapse_s.dta" // data set with 4950 cells, created in spatial_2_*  


scalar p = 0 // if 0, straight distance without decay : x^0 = 1, so all are equal.
scalar sigma = 47*0.59 // sigma at comp. line 
scalar kappa = 51.7*0.59 // kappa at comp. line 
scalar D = 3 // width of increment 

***************************
**** reduce dimensions ****
***************************

keep if year == 2011 & ccode == 1 // reduce data in memory 
  
sort cell // now sorted: first bunching (1-50), then restricted (51-1015) then rest. 

gsort -bunching +cell 
gen bcell = _n if bunching == 1 // identify B

gsort -restricted +cell 
gen rcell = _n if restricted == 1 // identify R


**********************************************
*** define upper and lower bounds of cells ***
**********************************************
		
	** gen mid-point for E 
gen a_mid = a_bin 
gen e_mid = `=kappa' + `=sigma'*a_mid //  

	** gen upper and lower level of E at mid-point (e.g. 6kg, 6.5kg ...)
replace e_bin = k_bin + a_bin*`=sigma' // gravity center: upper edge of mid-point 
gen e_ml = k_bin + a_bin*`=sigma' - `=D' // lower edge of mid-point 


replace eei_bin = (e_bin/(51.7+47*a_bin))* 100 // assert that bins are consistent with eei boundaries. 

	** label as safety net 
label variable e_bin "E_i at midpoint"

label variable a_min "a_i left"
label variable a_max "a_i right"
label variable k_bin "kappa (std. E)"
label variable e_mid "E_m"


keep eei_bin a_* e_* k_bin s_bin eei_bin cell bcell rcell // cut 

sort cell // order is crucial 

***********************************
******* start loops ***************
***********************************

**********************************	 
*  3a: neighborhood slope / fan *
**********************************

***********************************
*** gradients from unit circle  ***
***********************************
	
	** sw 
scalar g0 = tan((225*(_pi/180)))
scalar g1 = tan((202.5*(_pi/180)))
	** horizontal 
scalar g2 = tan((180*(_pi/180)))
	** mixed (nw) 
scalar g3 = tan((157.5*(_pi/180)))
scalar g4 = tan((135*(_pi/180)))
scalar g5 = tan((112.5*(_pi/180)))
	** vertical 
scalar g6 = -tan((90*(_pi/180))) // has to be minus due to quadrant switch at 90. 'approach' from left on unit circle.
scalar g7 = tan((67.5*(_pi/180)))
scalar g8 = tan((45*(_pi/180)))

forvalues i = 0/8 {
	scalar delta`i'= (1/g`i')
}	

* scalar list 

**************************************
* 3b: define scalar for every cell ***
**************************************


sum cell, meanonly 
scalar cmax=`r(max)'
local c = `r(max)' // scalar for each cell.  
	
sum bcell, meanonly 
scalar bmax=`r(max)'

sum rcell, meanonly 
scalar rmax=`r(max)'

	*** scale to normalize slope of compliance line to 1. From formula: EEI = (E / k + ba ) * 100
scalar sd_a = 100/`=sigma'
scalar sd_e = 100


forvalues i = 1/`=cmax' {
	* 5-digit cell-id
	scalar cell`i' = 10000 + `i'  
	
	* capacity 
	quietly sum a_bin if cell == `i'
	local a = r(mean)
	scalar a`i' = `a'/sd_a
	
	quietly sum a_min if cell == `i'
	local a = r(mean)
	scalar a_min`i' = `a'/sd_a
	
	quietly sum a_max if cell == `i'
	local a = r(mean)
	scalar a_max`i' = `a'/sd_a
	
	* std.  
	quietly sum s_bin if cell == `i'
	local s = r(mean)
	scalar s`i' = `s'
	
	quietly sum k_bin if cell == `i'
	local s = r(mean)
	scalar k`i' = `s'
	
	* eei 
	quietly sum eei_bin if cell == `i'
	local x = r(mean)
	scalar x`i' = `x'
	
	* e [four courners and midpoint]
	quietly sum e_bin if cell == `i'
	local e = r(mean)
	scalar e`i' = `e'/sd_e
	
	quietly sum e_low_left if cell == `i'
	local e = r(mean)
	scalar e_ll`i' = `e'/sd_e
	
	quietly sum e_up_right if cell == `i'
	local e = r(mean)
	scalar e_ur`i' = `e'/sd_e
	
	quietly sum e_up_left if cell == `i'
	local e = r(mean)
	scalar e_ul`i' = `e'/sd_e
	
	quietly sum e_low_right if cell == `i'
	local e = r(mean)
	scalar e_lr`i' = `e'/sd_e
	
	quietly sum e_mid if cell == `i'
	local e = r(mean)
	scalar e_mid`i' = `e'/sd_e
}


***********************************
**** 3c: define neighborhood ******	 
***********************************
		 
****************************	 
	* inclusion conditions 
	
forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k'
     local m = r(mean)
	  * m contains cellid of cell 1	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local n = r(mean)
		 * n contains cellid of cell 2 
		  forvalues d = 0(1)8 {
			** qualifying conditions: two parallel beams from edges of bunching cell touching the compliance line 
		  scalar qr`d'_`m'_`n' = (e_ur`m'+((a`n'-a`m')/delta`d')) // uppper-right boundary 
		  scalar ql`d'_`m'_`n' = (e_ul`m'+((a`n'-a`m')/delta`d')) // upperleft boundary 
          }
     }
}	 


******************************* 
***  4a: defining neighbors. ***
*******************************


forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k'
     local i = r(mean)
	  * i contains cellid of cell i in B 	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		  * j contains cellid of cell j in R 
		  
		  ** knock-out condition: use only relevant segment (moved up  for efficiency: saves 7*4950 scalars to compute)
		  scalar x_`i'_`j' 		= cond( x`j' > x`i'     			&  inrange(x`j', 59, 87)		,	1,0) // only neighbors in B. --> hadamard. 
		  
		  forvalues d = 0(1)8 {   
		   ** test boundaries against conditions : neighbor if gravity center falls in beam   
		   scalar s`d'_`i'_`j' 	= cond(ql`d'_`i'_`j' <= e`j' 		&  qr`d'_`i'_`j' >= e`j'		,	1,0)  // 
		   
          }
     }	
}


*********************************
*** 4b: transform to matrix *****
*********************************
display cmax // check dimension 

*** neighborhood matrices
mat define X  = I(`=cmax') // knock out criterium
mat define S1 = I(`=cmax') // sw 1   
mat define S2 = I(`=cmax') // sw 2
mat define S3 = I(`=cmax') // nw 1 
mat define S4 = I(`=cmax') // nw 2
mat define S5 = I(`=cmax') // nw 3
mat define S6 = I(`=cmax') // nw 4 
mat define S7 = I(`=cmax') // no 1
mat define S8 = I(`=cmax') // no 1

***************************
**  5: fill matrix ********
***************************

forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k' & bcell != . 
     local i = r(mean)
     forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l' 
          local j = r(mean)
		  ** run through combinations and put in distance **
		  mat X[`k',`l']=(x_`i'_`j')
		  mat S1[`k',`l']=(s1_`i'_`j')
		  mat S2[`k',`l']=(s2_`i'_`j')
		  mat S3[`k',`l']=(s3_`i'_`j')
		  mat S4[`k',`l']=(s4_`i'_`j')
		  mat S5[`k',`l']=(s5_`i'_`j')
		  mat S6[`k',`l']=(s6_`i'_`j')
		  mat S7[`k',`l']=(s7_`i'_`j')
		  mat S8[`k',`l']=(s8_`i'_`j')
         }
	 forvalues c = 1/`=cmax' {
	** replace own-cross with zero (identical cell)
	mat X[`c',`c']=0
	mat S1[`c',`c']=0
	mat S2[`c',`c']=0
	mat S3[`c',`c']=0
	mat S4[`c',`c']=0
	mat S5[`c',`c']=0
	mat S6[`c',`c']=0
	mat S7[`c',`c']=0
	mat S8[`c',`c']=0
        }
    } 

	
	** HADAMARD for each sector : eliminates those that are not from R or do not meet condition. 
forvalues i = 1/8 {
	matrix NS`i'=hadamard(X, S`i')
}	
	
***********************************************
****** 6: post results to weight matrix. ******
***********************************************
cd "$store"

	********************
	** Sectors 1-4 *****
	********************
preserve

svmat NS1, names(s1_)
svmat NS2, names(s2_)
svmat NS3, names(s3_)
svmat NS4, names(s4_)



save spatial_3_wmat_beam_1234, replace

restore


	********************
	** Sectors 5-8, Hadamard ******
	********************
preserve

svmat NS5, names(s5_)
svmat NS6, names(s6_)
svmat NS7, names(s7_)
svmat NS8, names(s8_)

	** check ** 
sum s5_82 s6_82 s7_82 s8_82 
tab cell if s5_82 > 0
tab cell if s6_82 > 0 
tab cell if s7_82 > 0 
tab cell if s8_82 > 0 

save spatial_3_wmat_beam_5678, replace

restore

log close 

clear 

exit 
	
	